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ABSTRACT 


This  report  presents  progress  on  an  ongoing  research  project  looking  at  near-field  infrasound  signals  as  a  basis  for 
discriminants  between  underground  nuclear  tests  (UGT)  and  earthquakes  (EQ).  In  an  earlier  program,  infrasound 
signals  from  UGTs  and  EQs  were  collected  at  ranges  of  a  few  hundred  kilometers,  in  the  far- field.  Analysis  of  these 
data  revealed  two  parameters  that  had  potential  for  discrimination  purposes:  signal  duration  and  wind-corrected 
amplitudes.  These  far-field  differences  should  be  present  in  the  near-field  signals  as  well.  To  study  the  near-field 
signals,  we  are  using  computational  techniques  based  on  modeled  ground  motions  from  UGTs  and  EQs.  One  is  the 
closed  form  numerical  integration  of  the  Rayleigh  integral  (RI),  and  the  other  is  the  application  of  a  time-domain, 
finite-difference  computational  fluid  dynamics  (CFD)  program,  CAVEAT.  This  report  summarizes  recent  progress 
in  modeling  these  signals,  showing  comparisons  of  waveforms  and  power  spectra  from  the  two  techniques.  There  is 
also  a  discussion  of  the  effects  of  spatial  and  temporal  zoning  on  the  quality  of  the  results.  Application  of  Fourier 
techniques  to  the  basic  ground  models  is  introduced  as  an  analytic  path  to  the  radiation  patterns  of  the 
ground-motion  sources. 
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OBJECTIVES 


The  objective  of  this  research  is  to  find  differences  in  the  near- field  infrasound  signals  of  UGTs  and  earthquakes  that 
can  be  the  basis  for  establishing  discriminants  between  the  two  sources.  Such  differences  in  the  near-field  signal 
would  likely  survive  to  a  longer  range. 

RESEARCH  ACCOMPLISHED 


We  have  continued  the  numerical  studies  on  modeling  the  near- field  infrasound  signal  from  the  surface  ground 
motion  from  underground  nuclear  tests  (UGT)  and  earthquakes  (EQ).  The  computational  tools  include  numerical 
integration  of  the  closed-form  RI  and  the  use  of  the  time-domain,  finite-difference  code  CAVEAT.  CAVEAT  was 
only  mentioned  in  last  year’s  MRR  report  but  has  been  used  more  extensively  this  past  year. 

CAVEAT  Developments 

CAVEAT  is  an  established  computational  tool  and  has  been  applied  to  calculations  of  the  time  evolution  of 
atmospheric  nuclear  explosions,  including  hydrodynamic  outputs  and  optical  signature  outputs.  CAVEAT  is 
documented  in  Adessio  et  al.  (1992).  The  following  succinct  description  of  the  CAVEAT  code  is  directly  quoted 
from  the  abstract  of  the  Adessio  et  al.  (1992)  report: 

CAVEAT  uses  an  explicit  time-marching,  conservative  finite-volume  numerical  technique  in  which  all 
state  variables,  including  velocity,  are  cell  centered;  values  at  vertices  and  cell  faces  are  derived.  The 
technique  is  a  variation  of  the  Godunov  method  that  uses  an  approximate  Riemann  solver  and 
accommodates  arbitrary  equations  of  state.  Spatial  differencing  may  either  be  first  order  (constant 
across  the  cell)  or  second  order  (linear  variation  across  the  cell)  with  a  choice  of  limiters  of  the 
gradient  in  an  attempt  to  preserve  monotonicity.  The  formulation  is  spatially  two-dimensional  with 
options  for  Cartesian  and  curvilinear  geometries,  e.g.,  cylindrical  (r,z).  Discretization  is  achieved  with 
a  mesh  of  arbitrary  quadrilateral  cells  whose  vertices  can  move  with  time.  Arbitrary  mesh  motion  is 
supported  by  allowing  transport  of  material  between  cells  according  to  the  arbitrary 
Lagrangian-Eulerian  (ALE)  technique.  The  computation  is  performed  in  two  phases:  a  Lagrangian 
phase  and  a  remapping  phase  in  which  conserved  variables  are  transferred  from  the  Lagrangian  mesh 
to  an  arbitrary  specified  mesh.  The  dynamic  mesh  capability  generally  smoothes  distortions  in  the 
mesh  and  can  also  result  in  higher  resolution  around  features  of  interest,  such  as  a  shock  discontinuity. 

The  report  contains  results  for  test  cases  of  a  shock  tube,  spherical  blast  wave,  a  pure  advection  problem,  and  a 
shock  on  wedge. 

The  initial  calculations  with  CAVEAT  used  a  separate  ground  motion  acceleration  model  that  was  not  the  same  as 
that  used  in  the  RI  code.  We  have  worked  to  correct  this  by  using  bi-cubic  spline  interpolations  on  the  RI  code 
computed  velocity  data  as  a  function  of  range,  r,  and  time,  t.  For  a  given  event,  the  full  set  of  surface  velocities, 
v(r,t)  is  written  to  an  input  file  to  CAVEAT.  At  the  bottom  of  the  CAVEAT  mesh  (grid),  we  employ  the  specified 
velocity  boundary  condition  whose  values  at  a  specific  time  and  range  are  interpolated  from  the  v(r,t)  field  using  the 
bi-cubic  spline  algorithms  from  Press  et  al.  (1990).  In  this  manner  the  velocity  data  in  the  RI  code  and  in  CAVEAT 
can  be  made  nearly  the  same.  This  scheme  is  quite  general  and  efficient  and  applicable  to  other  ground  motion 
sources  as  well. 

For  the  CAVEAT  results  reported  here,  we  used  cylindrical  (r,z)  coordinates  with  500  mesh  cells  in  each  direction. 
The  standard  run  was  with  radial  and  vertical  steps  of  30  m.  The  finer  zone  calculation  was  performed  with  mesh 
sizes  of  15  m.  The  runs  were  made  with  an  ALE  coefficient  of  0.9.  The  ambient  atmosphere  was  from  a  Committee 
on  Space  Research  (COSPAR)  International  Reference  Atmosphere  (CIRA)  standard  atmosphere. 

We  first  show  pressure  contours  for  the  standard  and  fine  zone  calculation  (Figure  1),  where  the  input  velocity  data 
were  from  the  modeled  accelerations  for  the  Tortugas  event  in  hole  U3gg.  In  the  CAVEAT  results,  we  show  the 
signal  pressure  that  is  the  total  pressure  minus  the  ambient  pressure,  sometimes  referred  to  as  the  perturbation 
pressure. 
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Figure  1.  Pressure  contours  for  the  standard  CAVEAT  run  (left)  and  finer  zone  calculation  (right)  at  20 
seconds.  The  x  and  y  values  are  in  centimeters,  and  pressure  contours  are  in  dynes/cm2.  Each 
calculation  used  250,000  zones.  Pressure  contour  values  are  shown  in  the  upper  right. 

The  finer  zone  result  shows  more  structure  in  the  pressure  field.  From  the  contour  legends  we  see  that  the  finer  zone 
result  has  larger  (in  absolute  value)  maximum  and  minimum  pressures.  By  taking  radial  slices  through  the  mesh  at 
different  elevation  angles  we  can  make  additional  comparisons.  An  example  is  shown  in  Figure  2. 
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Figure  2.  Radial  slices  through  the  Caveat  mesh  for  the  standard  run  (left  2x2  set,  A)  and  for  the  finer  zone 
run  (right  2x2  set,  B).  The  slices  are  at  elevation  angles  of  18,  26,  33,  and  45  degrees.  The  pressure 
unit  “mbar”  stands  for  “micro  bar.” 

The  finer  zone  run  generally  has  higher  frequency  structure  and  larger  amplitudes,  compare  the  3 3 -degree  panels.  To 
examine  the  impact  of  smoothing,  we  did  a  five-point  running  average  on  the  pressure  values  on  the  45 -degree  fine 
zone  data  (B  set),  and  the  result  is  quite  similar  to  the  45-degree  standard  run  result  (A  set).  These  results  show  the 
improved  velocity  model  that  has  now  been  implemented  into  CAVEAT. 
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Rayleigh  Integral  (RI)  Code  Developments 

For  review,  the  RI  is  given  by  the  following  equation: 


a(r,(t-R/ c))dA  ? 
R 


(1) 


where  a  is  the  acceleration  of  the  ground,  r  is  a  location  on  the  ground  (referenced  to  a  center  position),  dAs  is  an 
element  of  area  on  the  ground,  R  is  the  distance  from  the  ground  element  to  the  observing  point,  t  is  time,  p  is  the  air 
density,  c  is  the  speed  of  sound  in  air,  p  is  the  air  pressure,  and  Rq  is  the  observation  location. 


The  RI  code  uses  modeled  ground  motions  for  some  30  events.  The  parameters  include  the  peak  acceleration  and 
temporal  width  of  the  initial  acceleration  pulse  and  the  same  for  additional  contributions  from  plastic  or  cavity 
rebound  signals.  A  specified  exponential  decay  with  range  is  also  applied. 


Figure  3.  An  illustration  of  the  characteristics  of 
ground  motion  from  a  UGT. 


Figure  3  illustrates  the  features  in  UGT 
acceleration  records.  Some  events  have  just  the 
initial  and  slapdown  contributions,  such  as 
Tortugas,  U3gg,  while  others  have  additional 
contributions  such  as  Duoro,  U31v.  When  the 
acceleration  is  the  simple  two-pulse  shape,  the 
velocity  exhibits  a  simple  form  as  well.  This  is 
illustrated  in  Figure  4.  The  acceleration  is  a  classic 
two-pulse  record,  and  the  velocity  has  a 
well-defined  linear  portion  due  to  the  -1  g  spall 
phase  in  the  acceleration,  Jones  et  al.  (1993)  and 
Huan  and  Walker  (1980). 
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Figure  4.  A  measured  ground-motion  acceleration  record  for  Totugas  (left)  and  the  velocity  record  (right), 
where  units  are  gs  for  acceleration  and  m/s  for  velocity.  The  horizontal  axis  is  the  sample  number. 

If  the  acceleration  record  has  contributions  from  other  components,  the  behavior  is  more  complicated,  as  is 
illustrated  in  Figure  5  for  the  Duoro  event. 
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Figure  5.  The  same  as  Figure  4  except  for  the  Duoro  event,  which  shows  more-complicated  behavior. 

We  use  the  modeled  ground  motions  from  these  two  events  to  illustrate  some  results  from  the  RI  code.  First,  we 
show  computed  power  spectra  for  the  RI  code  calculation  for  an  observer  at  an  horizontal  range  of  5  km  and  a  height 
of  4  km  and  for  an  observer  at  a  horizontal  range  of  15  km  and  a  height  of  10  km.  These  are  close  to  the  same 
elevation  angle.  Figure  6  presents  this  comparison  and  shows  good  agreement. 


PSD  U3ggz4f5  PSD  u3gg  zl  Orl  5 


Figure  6.  Power  spectra  (relative  units)  for  the  Tortugas  event,  hole  U3gg,  for  a  height  of  4km  and  a  range  of 
5  km  (left)  and  for  a  height  of  10  km  and  range  of  15  km  (right).  The  results  are  reassuringly  close. 
Frequency  is  on  the  horizontal  axis. 

To  illustrate  the  effects  of  the  more  complicated  ground  motion,  Figure  7  shows  the  RI  code  power  spectra  for  the 
Duoro  event.  As  expected,  the  Duoro  event  has  larger  contributions  at  higher  frequencies  than  does  the  Tortugas 
event,  due  to  the  more-complicated  source  acceleration. 
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PSD  U3tVZ4*5 


Figure  7.  RI  code  power  spectra  (relative  units)  for  the 
Duoro  event  at  a  range  of  5  km  and  a  height 
of  4  km.  Frequency  is  on  the  horizontal  axis. 
Compare  this  result  with  that  in  the  left  panel 
of  Figure  6. 


Initial  Comparisons  of  Two  Techniques 

Below  we  show  a  comparison  between  the  CAVEAT  and  RI  codes  for  the  Tortugas  event,  U3gg.  CAVEAT 
advances  in  time  so  that,  at  a  given  time,  we  can  take  a  snapshot  of  the  pressure  field  in  the  computational  domain. 
Then  we  can  plot  radial  slices  through  that  domain  at  different  elevation  angles.  This  is  shown  in  Figure  8  with  the 
left  four  plots  from  CAVEAT.  The  RI  code  provides  a  time  history  at  a  given  observation  location.  When  looking  at 
the  CAVEAT  radial  slices,  one  must  remember  that  the  pulse  is  proceeding  to  the  right.  Thus  to  properly  compare 
the  two  code  results  one  should  flip  the  CAVEAT  plot.  With  that  in  mind,  one  can  see  that  the  waveforms  are  quite 
similar  in  shape.  The  RI  code  result  is  at  45  degrees  and  shows  a  little  better  agreement  with  the  CAVEAT  profile  at 
26  or  18  degrees.  This  is  probably  due  to  refraction  in  a  layered  atmosphere  as  is  used  in  CAVEAT  while  a  uniform 
atmosphere  is  used  in  the  RI  code. 
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Figure  8.  A  comparison  of  CAVEAT  and  RI  code  results  for  the  Tortugas,  U3gg,  event.  The  left  four  panels 
are  CAVEAT  profiles  at  20  seconds.  The  CAVEAT  results  are  for  the  finer  zoned  calculation.  The 
right  plot  is  the  RI  code  result  at  45  degrees  and  range  of  6.5  km.  Due  to  a  plotting  error,  the  RI 
code  pressures  should  be  a  factor  of  10  larger. 
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Some  Analytic  Considerations 


1.  An  obvious  direction  in  this  work  is  that  of  a  Fourier  approach  to  the  RI  formulation.  This  is  would  serve  as  an 
easy  way  to  compute  the  radiation  pattern  of  ground-motion  sources.  A  nice  approach  was  given  by  Wecksung 
(1985)  in  an  unpublished  report  and  is  summarized  below  for  ease  of  reference,  with  some  notational  changes. 

The  basic  integral  is  given  by 

2^  A  r 

The  area  of  motion  is  smaller  than  R0  the  distance  from  the  center  of  the  motion  to  the  observation  point.  The 
source  point  on  the  ground  is  at  ro  =  ro(x,y,0),  and  r  is  the  distance  from  the  source  point  to  the  observation  point  R. 
It  is  possible  to  use  1/ R0  for  the  1/r  term  in  the  integral  and  take  it  outside.  Apply  a  change  of  variable  of  r  =  t-r/c 
and  take  a  temporal  Fourier  transform  of  Eq.  (2)  to  get 


p(R9G))- 


Po 

2nR, 


Jj*  J  ( a(x,y,r)dxdy)e  l27rTdre~ 


ilnrlc 


0  A  -oo 

The  normal  acceleration,  a,  is  related  to  the  normal  velocity,  u,  by 

du(x,y,t) 


a{pc,y,t)  —  ■ 


dt 


and  the  integral  over  a(x,y,t)  becomes 

00 

|  a(x,y,  z)e~lcox dr  =  icou(x,y,co)^ 

-00 


(3) 

(4) 

(5) 


where  we  have  replaced  27if  with  co  and  u  is  the  temporal  Fourier  transform  of  u.  Then  using  Eq.  (5),  we  can  write 

r 

p(R,co )  =  —  \\ u(x,y,a>)e~,wr/cdxdy  ■  (6) 

IttR, 

At  this  stage  the  basic  acceleration  integral  has  been  transformed,  through  a  temporal  Fourier  transform,  to  a  spatial 
Fourier  transform  of  the  normal  velocity  transform.  Next,  let  e  be  a  unit  vector  in  the  direction  of  R, 
e  =  (sin  #  cos  sin  #  sin  cos  0) .  It  is  possible  to  write  r  =  R-ro  and  from  that,  have  r2  =  R2  -  2R  •  r  +  r02 .  A  little 
simplification  gives  an  approximate  expression  for  r,  r  «  R-  (xsin#cos^  +  ysin#sin^) .  Next,  for  the  spatial 
transform  of  u  ,  write 


Ua(vx,vy)  =  JJ  u(x,y,6))exp[-i27r(xvx  +yuy)\dxdy  ■ 

s 

and  with  the  approximate  expression  for  r,  Eq.  (6)  may  be  rewritten  as 

-^sin^cos^  -co sin# sin 


\  iO)P0  /  iCDRXTT 

p(R,  a)  =  y—x  exp( - )Ua 

2nR  c 


2  nc 


2  nc 


y 


The  acoustic  intensity  is  related  to  the  pressure  by 
r  \p{R,cof 

Ia(R)  =  - - L’ 

P0c 

which  leads  to  an  expression  for  the  far  field  intensity 


T  fL_  (ML 

J)  4  tt2R2c 


U„ 


-fflsin^cos^  -fflsin^sin^  j 


2  nc 


Inc 


(7) 

(8) 


(9) 


(10) 
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2.  Below  is  given  a  simple  and  straightforward  derivation  of  an  expression  for  the  maximum  amplitude  of  a  uniform 
piston.  Let  the  piston  displacement  be  given  by 

z  =  Asin(27rft) ,  (11) 

where  A  is  the  amplitude  of  the  motion,  f  is  frequency,  and  t  is  time.  Then  the  velocity  is  given  by 

V  =  27ifAcos(27ift) .  (12) 

On  axis,  the  solution  for  pressure  from  an  uniform  piston  is  proportional  (Pierce  1989)  to  a  difference  of  velocities  at 
retarded  times, 


p(h)  =  pc 


c  c 


=  pc  A  V  ’ 


(13) 


where  h  is  the  height  above  the  center  of  the  piston,  R  is  the  radius  of  the  piston,  p  is  the  air  density,  and  c  is  the 
speed  of  sound.  From  Eq.  (12),  the  AV  term  can  be  written  as 


A  V  =  Aa 


h  \[h 2  +i?2 

cos (a(t - ))  -  cos (a(t  - - - )) 

c  c 


(14) 


where  a  =  27if.  We  assume  h  »  R,  then  the  square  root  term  becomes 


h  R2  h 

-  +  —  =  —  +  y 

c  2ch  c 


(15) 


Let  (t-(h/c))  =  b,  substitute  Eq.  15  in  Eq.  14,  expand  the  second  cosine  term,  cancel  obvious  terms  and  use  the  small 
angle  approximations  for  sine  and  cosine  to  get 


AV  =  -Aa2ysin(ab) , 

with  a  maximum  of  Aay.  Then  Eq.  (13)  can  be  rewritten  as  (for  the  maximum  value) 


p{h)  =  pcAa2y  = 


2  n2Apf2R2 
h 


(16) 


(17) 


Evaluation  then  requires  only  the  displacement  amplitude,  frequency,  piston  radius,  air  density  and  altitude.  An 
equivalent  expression  to  Eq.  (17)  was  given,  without  derivation,  in  Banister  (1979)  but  had  been  overlooked  by  the 
author. 


Some  Additional  Items 

Ground-motion  data,  from  measured  ground-motion  data  files,  for  the  30  modeled  events  have  been  extracted.  The 
accelerations  and  velocities  plots  have  been  put  into  a  draft  report  that  can  be  supplied  upon  request.  Power  spectra 
are  now  routinely  computed  for  any  case  calculated  by  the  RI  code.  A  postprocessor  was  written  to  extract  radial 
slices  from  the  CAVEAT  mesh. 


CONCLUSIONS  AND  RECOMMENDATIONS 


The  improved  velocity  model  from  the  RI  code  input  to  CAVEAT  was  implemented  just  before  the  deadline  for  the 
submission  of  this  paper;  only  a  few  runs  have  been  made.  The  mesh  resolution  in  the  CAVEAT  runs  needs  further 
exploration.  Initial  comparisons  of  CAVEAT  profiles  with  the  RI  code  show  some  differences  that  need  to  be 
studied  and  understood. 
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